
n = 300;
X = randn(n,10,'single');

effect = [zeros(5,1); 0.5*(-1).^(rand(5,1)<.5)];
effect_prec = [randn(1,1); zeros(4,1); randn(2,1); zeros(3,1)];

mu = 1./(1+exp(-X*effect-2));
p = exp(X*effect_prec+1);

a = mu .* p;
b = (1-mu) .* p;

y = betarnd(a,b);

[slab0,scale_delta0,llik0] = sample_beta_reg(y,X,100,1000);

[spike,slab,prec_spike,prec_slab,llik] = sample_spike_slab_beta_reg(y,X,repmat(0,size(X,2),1),100,1000);



figure(1);
subplot(2,1,1);
boxplot((spike.*slab)'); hold on; plot(1:numel(effect), effect, 'gd', [0,(1+numel(effect))], [0,0,], 'k--'); hold off;
subplot(2,1,2);
boxplot(slab0'); hold on; plot(1:numel(effect), effect, 'gd', [0,(1+numel(effect))], [0,0,], 'k--'); hold off;


figure(2); boxplot((prec_spike.*prec_slab)'); hold on; plot([0,(1+numel(effect))], [0,0,], 'k--'); hold off;
figure(3); hist(y,20);


[r,p] = corr(X,y,'type','Spearman');
[~,p_bet] = ttest(slab');
[~,p_bet0] = ttest(slab0');

p(p < 1e-10) = 1e-10;
p_bet(p_bet < 1e-10) = 1e-10;
p_bet0(p_bet0 < 1e-10) = 1e-10;

figure(4); plot(1:numel(p), -log10(p), 'ko', ...
                1:numel(p_bet), -log10(p_bet), 'r+',...
                1:numel(p_bet0), -log10(p_bet0), 'b^',...
                'linewidth', 2);

figure(5); semilogx(llik)

